{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Chapter 4 Exercises\n",
    "\n",
    "In this notebook we consider a subset of problems to do some calculations and shed light as to how do certain proofs / exercises"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "from numpy.random import multivariate_normal, seed\n",
    "from numpy.linalg import inv\n",
    "from sklearn.preprocessing import OneHotEncoder"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "%config InlineBackend.figure_format = \"retina\"\n",
    "np.set_printoptions(precision=4, suppress=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4.2\n",
    "Consideer the minimization of a sum-of-squares error function given by\n",
    "\n",
    "$$\n",
    "    E_D\\left(\\tilde{\\bf W}\\right) = \\text{Tr}\\left([\\tilde{\\bf X}\\tilde{\\bf W} -  {\\bf T}]^T[\\tilde{\\bf X}\\tilde{\\bf W} -  {\\bf T}]\\right)\n",
    "$$\n",
    "\n",
    "and suppose that all of the target vectors in the training $\\{{\\bf t}_n\\}_n$ set satisfy a linear constraint\n",
    "\n",
    "$$\n",
    "    {\\bf a}^T{\\bf t}_n + b = 0\n",
    "$$\n",
    "\n",
    "Show that as a consecuence of this constraint, the elements of the model prediction ${\\bf y}({\\bf x})$ given by the least squares solution\n",
    "\n",
    "$$\n",
    "    \\tilde{\\bf W} = \\left(\\tilde{\\bf X}^T\\tilde{\\bf X}\\right)^{-1}\\tilde{\\bf X}^T{\\bf T}\n",
    "$$\n",
    "\n",
    "also satisfies this constraint so that \n",
    "\n",
    "$$\n",
    "    {\\bf a}^T{\\bf y}({\\bf x}) + b = 0\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# According to Bishop, a clear example of this is by one-hot-encoding all target variables\n",
    "# so as to satisfy a linear constraint of the form 1^t x + 0 = 1 (p.145)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.axes._subplots.AxesSubplot at 0x1a20d91be0>"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAxcAAAIPCAYAAAAB0iN+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAWJQAAFiUBSVIk8AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzs3Xd8leX9//HXdZ/sCSHsFQgbBQTZLqDiat2jWneXdeDqzw7biv3Wftt+66JqrbV1tmod1VoVsSAOpoSCisxAgMgIm5CQkJz7+v1xcmISMg7JyZnv5+ORx4H7/pz7XIcAud/nWsZai4iIiIiISFs54W6AiIiIiIjEBoULEREREREJCoULEREREREJCoULEREREREJCoULEREREREJCoULEREREREJCoULEREREREJCoULEREREREJCoULEREREREJCoULEREREREJCoULEREREREJCoULEREREREJCoULEREREREJCoULEREREREJCoULEREREREJCoULEREREREJioRwN0CaZozZBGQBRWFuioiIiIjEtjzgoLW2X1suonAR2bJSU1Nzhg4dmhPuhoiIiIhI7Fq9ejWHDx9u83UULiJb0dChQ3MKCgrC3Q4RERERiWFjxoxh+fLlRW29juZciIiIiIhIUChciIiIiIhIUChciIiIiIhIUChciIiIiIhIUChciIiIiIhIUChciIiIiIhIUChciIiIiIhIUGifixjhui579+6ltLSUyspKrLXhblJcM8aQnJxMZmYmOTk5OI5yvIiIiMQ+hYsY4LouW7dupby8PNxNkRrWWioqKqioqKCsrIzevXsrYIiIiEjMU7iIAXv37qW8vJyEhAS6detGenq6bmTDzHVdysrK2LFjB+Xl5ezdu5fc3NxwN0tERESkXekONAaUlpYC0K1bNzIzMxUsIoDjOGRmZtKtWzfgq++RiIiISCzTXWgMqKysBCA9PT3MLZGG/N8T//dIRERih+vWn9/oar6jiIZFxQL/5G31WEQeYwyAJtiLiMQQ17U4jqF4XzmvLC9mz6EjdMlM5qIxvejVMa32vEg8UrgQaUf+cCEiIrHB61qqvS53vriStz7bTt3Pjh6au57zRvbkdxePIAEUMCQuKVyIiIiIBMjjGL7zzHLeX1ty1Dlr4fUVX1JZ7eWPV44JQ+tEwk/jaEREREQC4HUtyzfvazRY1PXO5zv4YtsBvK6GxEr8UbgQERERCYDHMbzwyZaAav++dCseDYuSOKRwITGlqKgIYwzXXnttuJsiIiIxqHjv4cDq9mljW4lPChciIiIiAcpKDWy6alZKYju3RCQyKVyIiIiIBMDrWr4+okdAtd8Y2V1zLiQuKVxIzJg5cyb9+vUD4JlnnsEYU/v19NNPh7dxIiIS9TyO4azju9EnJ63ZuvzOGUwb0lVzLiQuaSlaiRmnnXYa+/fv5+GHH2bkyJGcf/75tedGjRoVxpaJiEisMBievX4c33pyCV/uP3r+RZ+cNJ65fizqs5B4pXAhMeO0004jLy+Phx9+mFGjRjFz5sxwN0lERGKMxzH0zknjvTtO4bXlX/La8mJ2HzpC58xkLhrdiwtH9yTR46jXQuKWwoWIiIjIMfA4hrSkBK4Y34crJ/StPe66VrtyS9zTnAsRERGRVnBM/SChYCGicCEiIiIiIkGicCEiIhLnXFt/+rGrJVRFpJU050JiisfjAcDr9Ya5JSIikc+1FscYdpVW8kpBMTsOVNAxPZELTuhFv9x0vK7VxGQROSYKFxJTOnbsiDGGLVu2hLspIiIRzXUtLpaf/fNzXlq2td6Gb7PmbuCM4V158LJRJCd4FDBEJGAKFxJTMjIyGD9+PB999BHf+ta3GDRoEB6Ph3PPPZcRI0aEu3kiIhHDcQw/euVTXl5W3Oj5d1ft5NvPLOOF704IcctEJJrF3ZwLY0wnY8x3jDH/NMZsMMYcNsYcMMZ8bIz5tjHmmP5MjDG9jDF/NcZsM8ZUGmOKjDEPGWM6ttd7kOY999xznHPOOcyePZt7772Xn//85yxfvjzczRIRiRiuaynaXcYrBY0HC79FhXv4eP2uer0aIiLNiceei0uAPwLbgfeBLUBX4ELgSeAsY8wl1toW/yc1xuQDC4EuwBvAGmAccCtwpjFmsrV2T7u8C2nSgAEDePPNN8PdDBGRiOU4hheWbqHln3TwwtKtnDSwc/s3SkRiQjyGi3XAucBb1lrXf9AY81NgKXARvqDxagDXegxfsJhhrf1DnWs9ANwO3AfcELymi4iIBEfxvsMB1W3dV97OLRGRWBJ3w6KstfOstW/WDRY1x3cAj9f89rSWrmOM6Q9MB4qARxucvgcoA64yxqS3tc0iIiLBlpkS2OeLgdaJiEAchosWVNU8VgdQO7XmcU4jQaUUWACkAZoJJyIiEcXrWr4+okdAtV8f0SNq9r3wuhava6l2XQIY3Swi7UAfR9QwxiQAV9f8dnYATxlc87iuifPr8fVsDALmtvDaBU2cGhJAO0RERI6JxzGcNDCX4T2yWLXtYJN1XTKTufCEnjgRvhSttRZrYf66ErbsKSc9OYHpw7rSIS0J17UR336RWKJw8ZXfAMcBb1tr3w2gPrvm8UAT5/3HO7S1YSIiIsHmdS1/vXYslz+xmI27y4463zkjmee+PY4ET+QPcvjnf7/k/95dy/YDFbXHfp7gcP4JPfnlucNJwNFeHSIhonABGGNmAHfiW+3pqmBdtuaxxX5Za+2YJtpVAIwOUntERERqeRxDbkYy79x6Mv9c8SUvLytm58EKOqQlcv6onlw2tjdpSQkRf1P+0idb+dGrnx51vLLa5aVPtlK8r5xnrx8fhpaJxKe4DxfGmJuAh4EvgGnW2r0BPtXfM5HdxPmsBnUiIiIRxeMYPI6HS8f05ptj+9Qed63FMZEdKgAqq738+u3VzdYs2LCHd1ft4Izh3SI+KInEgsjv62xHxpjbgEeAz4EpNStGBWptzeOgJs4PrHlsak6GiIhIRGg4JyEagoXXtfx75XYOHK5qsfZvSzYrWIiESNyGC2PMj4AHgRX4gkXJMV7i/ZrH6Q139TbGZAKTgcPA4ra2VUREROrzOIZ1O0sDql2341A7t0ZE/OIyXBhjfo5vAncBvqFQu5upTTTGDKnZjbuWtbYQmAPkATc1eNq9QDrwrLX26FlyIiIi0mbJiYHdxgRaJyJtF3dzLowx1wC/BLzAR8AMc3T3b5G19umaX/cEVgOb8QWJum4EFgKzjDHTaurGA1PwDYe6O/jvQERERLyu5WtDuzJr7oYWa6cN7RqCFokIxGG4APrVPHqA25qo+QB4uqULWWsLjTEn4gsrZwJnA9uBWcC9xzA5XERERI6BxzGM6NWB0X06sHzL/mbrrpnYV/tdiIRI3PUTWmtnWmtNC1+n1akvqjmW18T1tlprr7PWdrfWJllr+1prb1WwEBERaV9e1/LIFaPp1TG10fMex/Dbi46nf+cMBQuREInHngsRERGJAR7H0DUrhbdmnMyzi4p4celWvtx/mOQEh7OP7871k/M4vpf2shUJpbjruZDYV1xczPXXX0+PHj1ITk4mLy+P2267jX379h3Tdfbu3cttt91GXl4eycnJ9OjRg+uvv57i4uJ2armIiBwrj2PISkng5ikDWPDjqWy47yzW/uosHrxsFMN6NLUVlYi0F/VcyDFZt7OUBRt2c6iimoyUBCYPyGVQ18xwN6tWYWEhkyZNoqSkhPPOO48hQ4awdOlSHn74YWbPns2CBQvo1KlTi9fZs2cPkyZNYt26dUydOpVvfvObrFmzhqeeeoq33nqLRYsW0b9//xC8IxERaUndhVkSPF99bqq9LURCT+FCArJgw24enruepZuOnkoyrl8Ot04byOQBuWFoWX033ngjJSUlzJo1i1tuuaX2+B133MGDDz7I3XffzeOPP97idX7605+ybt06br/9dh544IHa47NmzeLWW2/lxhtvZPbs2e3yHkREJD55XVsbiOr+WiSaGGttuNsgTTDGFIwePXp0QUFBs3WrV68GYOjQoe3Sjpc+2cJPXvsMt5m/Ko6B31w4gkvH9m6XNgRi48aN5Ofnk5eXR2FhIY7z1adXpaWldO/eHWstJSUlpKenN3mdsrIyOnfujOM4bN++nczMr3pmXNclPz+foqIiCgsLA+q9aO/vj4iIRDd/kFiz/SDvfrGTiiNeBnTN4OsjupOc4Al38yROjBkzhuXLly+31o5py3U050KatWDD7haDBYBr4cevfcqCDU3uR9ju5s2bB8D06dPrBQuAzMxMJk+eTHl5OYsXN79p+qJFizh8+DCTJ0+uFywAHMdh+vTpALz//vuNPV1ERCRgXteyt6ySy/60iDMf/ogH31vHHz8o5M5/rGTcfXP5+5It4W6iyDFRuJBmPTx3fYvBws+1MGvu+vZtUDPWrl0LwKBBgxo9P3DgQADWrVsXkuuIiIg0x1pLRZWXS/+0mCWNDDs+cLiKn/7zM/6xbGsYWifSOgoX0qR1O0sbnWPRnCWb9rJuZ2k7tah5Bw4cACA7u/HVQfzH9+9verOlYF5HROKLt8EnMdZaNPRYWvL84s1s2l3WbM3v311LteuGqEUibaMJ3dKk1g5xWrBhd0StIOXn/yFfd1WRcF5HRGKDay2OMazefpDnFm9mzfaDJCU4TBnchSvG9yE7NVH/X0ijjDG8+EnLvRIlpZXMXV3CtKFdSHD0ubBENoULadKhiuqQPq+t/D0K/p6Hhg4ePFivrr2vIyKxz/9hw49e/ZSXGtwkflK0j1nz1vPI5aP52rCu4WieRIGiPc33Wvht2l2Go5AqUUDxV5qUkdK67Nna57XV4MGDgabnQqxf75sP0tRcimBfR0RinzGGh/6z7qhg4VdR5XLT35ezdkfpUcOmRADSkwL7mZmRnIBG2Uk0ULiQJrV234pw7XcxZcoUAObMmYPbYGxqaWkpCxYsIDU1lQkTJjR7nQkTJpCamsqCBQsoLa0/f8R1XebMmVPv9UQkflVWeXl6YVHzNdUuT360UXsWSKPOOr5bizUJjuHM47qhv0ISDRQupEmDumYyrl/OMT1nfL+csM23yM/PZ/r06RQVFfHoo4/WO3fPPfdQVlbG1VdfXW+PizVr1rBmzZp6tRkZGVx11VWUlZUxc+bMeuceeeQRioqKOOOMM7RDt0icq3Zd5q4p4eDhloeC/mvlNlz1XEgDXtfy7cn9SGghNXxjZA9yM5I1d0eiguZcSLNunTaQq/6yJKDlaB0DM6YNbP9GNeOxxx5j0qRJzJgxg7lz5zJ06FCWLFnC+++/z6BBg7jvvvvq1fs3tmu4osuvf/1r5s+fzwMPPMCKFSsYN24cq1ev5o033qBLly5HhRcRiU/7y48EVFdZ7VJR7SUtwCEwEh88jmFI9yx+d/EI7nrlU6ob+WE7Nq8jv77geO3YLVFDPRfSrMkDcvnfC49vsSvWv0N3uIZE+eXn57Ns2TKuvfZalixZwv33309hYSEzZsxg0aJFdOrUKaDrdOrUiUWLFjFjxgw2bNjA/fffz5IlS7juuusoKCggPz+/nd+JiEQ6A3TNTgmoNislgZRE7bQsjbtwdC/m3H4K10zKo0d2CjnpSYzrl8NDl43ihe9OICnBUbCQqGG0BnfkMsYUjB49enRBQUGzdatXrwa++hS+PSzYsJtZc9c3usnP+H45zJg2MOzBIlKF4vsjIuFR7bpM+t95lJRWNlt33eQ87vnG8BC1SqKR61qcBgFCvRUSSmPGjGH58uXLrbVj2nId9c9KQCYPyGXygFzW7SxlwYbdHKqoJiMlgckDciNyTwsRkVBIcBxumTqQn7/xeZM1WakJfOek/o3ePIr4NfZ3Q8FCopHChRyTQV0zFSZEROq4amJfDlZU8dB/1lHlrT8aoFtWCk9cPYbuHVK0R4GIxAWFCxERkTa6acoALh/XhxeWbmHtjlISPYapQ7pwxnHdcDAKFgHwuhZrLcaAR7tQi0QthQsREZEg6JiWyI2n5dcuF6rx8oGx1uJa+HjDLkoOVpKbkczJg3JJUMAQiUoKFyIiIkHQcA8CBYvAPLNoM3+cv4GdB7+aFN85M5kbTu3Pt0/qX9OboT9LkWihcCEiIiJh8eB763h47vqjju8qreR//r2a3aVH+NFZQ8LQMhFpLfU5ioiISEi51rL9wGH+MO/oYFHX4x8WsnVvOa6WzReJGgoXIiIiElKOMfxtyRYa2ZC6Hmvh2UWbNSFeJIooXIiIiEjIrfryYEB1X2w/0M4tEZFgUrgQERGRkEtKCKw3IlGrRolEFf2LFRERkZByXcupgzoHVHvq4M6acyESRRQuREREJKQcx3Dh6F5kpyY2W5eRnMClJ/bWnAuRKKJwISIiIiGXlODw2LdGk5roafR8coLDo1ecQFpS4+dFJDIpXIiIiEjIOcYwKb8Tb94ymUtP7E1Kou+WJDnB4eIxvXjzlpM4ZVBnbaAnEmUULiSmvPLKK9xyyy2cfPLJZGVlYYzhyiuvbNW1iouLuf766+nRowfJycnk5eVx2223sW/fviC3WkQkPhlj6Jebwe8uHsHnM89gxS9O5/N7z+D3l4wkv3OGgoVIFNIO3XJsSlbDxg+gshSSM6H/qdBlaLhbVetXv/oVK1euJCMjg169erFmzZpWXaewsJBJkyZRUlLCeeedx5AhQ1i6dCkPP/wws2fPZsGCBXTq1CnIrRcRiT8exxcgEjwOHdKSjjouItFF4UICs3E+fPA72Lzg6HN9J8Opd0H/00LcqKM9+OCD9OrViwEDBvDBBx8wZcqUVl3nxhtvpKSkhFmzZnHLLbfUHr/jjjt48MEHufvuu3n88ceD1WwRERFpgWvtUZP7va5VEI0wGhYlLVv+LDx3QePBAnzHn7sAlj8X2nY1YsqUKQwcOLBNXekbN25kzpw55OXlcdNNN9U7d++995Kens5zzz1HWVlZW5srIiIiAbDWsn1/Bb9+ezUn/XYeJ/xyDhc8uoBXCoo5Uu3ibWm7dwmZuAsXxpiLjTF/MMZ8ZIw5aIyxxpjnW3GdoprnNva1oz3aHhYb58Obt4J1m6+zLrw5w1cf5ebNmwfA9OnTcRps3pSZmcnkyZMpLy9n8eLF4WieiIhIXLHWMm9NCVPvn88TH26keN9h9pVX8d+t+/nRq59ywWMLOFRRpYARIeIuXAA/A24GRgFftvFaB4B7G/n6fRuvGzk++F3LwcLPuvDB/7Vve0Jg7dq1AAwaNKjR8wMHDgRg3bp1IWuTiIhIPHKtZduBCm7823Iqqxu/H1m17SAzXlyh4VERIh7nXNwOFAMbgFOB99twrf3W2pnBaFREKlnd9FCopmz+2Pe8CJrkfawOHDgAQHZ2dqPn/cf3798fsjaJiIjEI8cYnl1Y1GSw8Ptg3S4Kdx2iX6d0HIWMsIq7ngtr7fvW2vXWWvWdtWTjB6F9XpTw/9XREokiIu3DWou1FtdaNu8pY+OuQxypubl09eM77rzzeWCjzf/96XYFiwgQjz0XwZRsjLkS6AOUAZ8CH1prveFtVpBUlob2eRHC3zPh78Fo6ODBg/XqREQkeLyuxWJ58sNNPL9kM8X7DgPQIS2RS8b05papA0hPTtAQmDhSVlkdUF15gHXSvhQu2qYb0HCJpE3GmOustQF/fG+MKWji1JBWtywYkjND+7wIMXjwYKDpORXr168Hmp6TISIirWOtL1h8++llfLBuV71z+8ur+PNHG5m3poTXfjCRjJREBYw40bdTOnvKjgRU57pWvRdhFnfDooLoKWAavoCRDhwP/AnIA94xxowMX9OCpP+poX1ehPDvjTFnzhxct/4Yz9LSUhYsWEBqaioTJkwIR/NERGKWBZ5eUHRUsKircNch/uet1QoWccLrWi4f17vFuvQkD+ef0EPBIgIoXLSStfZea+08a+1Oa225tfZza+0NwANAKjDzGK41prEvoHXbSwdLl6G+DfKORd+TomYyd1VVFWvWrKGwsLDe8fz8fKZPn05RURGPPvpovXP33HMPZWVlXH311aSnp4eyuSIiMc8xhucWb26x7s2V2zhwuCoELZJw8ziG80b15LieWc3W3TJtIGlJGpATCfRdCL7HgTuBU8LdkKA49S7fBnmBLEdrHDj1/7V/m5rx+uuv8/rrrwOwY4dvAtiiRYu49tprAcjNzeX3v/etFPzll18ydOhQ+vbtS1FRUb3rPPbYY0yaNIkZM2Ywd+5chg4dypIlS3j//fcZNGgQ9913X8jek4hIvNhxoILNe8pbrKusdvmkaC9TBndRD0Yc8DiGv31nAre9uIL315bUO5eRnMDNUwdww6n5YWqdNKRwEXz+v/Wx8bF2/9PgGw+3vJGeceAbs3z1YbRixQqeeeaZesc2btzIxo0bAejbt29tuGhOfn4+y5Yt4xe/+AWzZ8/m7bffpnv37syYMYN77rmHnJycdmm/iEg8swS+ElQ0LfroX/kKjMJQK3gcQ0ZyAk9dN5ZNu8v496fbKK/00qdTGueN6qEeiwij70bwTax53BjWVgTT6KuhQx/fBnmbPz76fN+TfD0W/U8LdcuOMnPmTGbOnBlQbV5eXrM/nHr37s1TTz0VpJaJiEhLumSm0DUrmZ0HK5utS3AMo3p3INLv072uxeMYivaUsXbHIZITHCbmdyIl0aOJx8fIH8r65qRxy1TfZrautThaFj7iKFw0wxiTCOQDVdbawjrHhwPbrbV7G9T3BR6p+e3zIWtoKPQ/zfdVstq3j0VlqW9VqP6nRs0cCxERiWyOgcvH9eGh/6xvtu5rw7rSOTMlRK1qHWstq7Yd4Ndvr2bxxq9uFzKTE7h0bG/uOnMwCTjqyThGdQOZgkVkirtwYYw5Hzi/5rfdah4nGmOervn1bmvtD2t+3RNYDWzGtwqU3yXAj40x7wObgFJ8IeQcIAV4G2h57E006jJUYUJERNrN90/JZ+7qEj77svG9hrpkJvOLrw+r7RWIRK5r+bR4P9/882IqqhqsOlhZzV8+3sTq7Qd59vpxQGS+B5HWirtwAYwCrmlwrH/NF/iCxA9p3vvAYOAEfMOg0oH9wMf49r14TjuAi4iIHBtjDEkJDi9+bwK/nb2G15Z/yaGajdESPYYzh3fjJ2cPpVtWSkQPKXIcw0//+flRwaKuhYV7eGV5MZeO6R3R70XkWMVduLDWziTAZWKttUU08pFCzQZ5AW+SJyIiIoHxOIbURA+/PO84fnzWEP67ZT9e13Jcz2xy0pPwRvhcBa9rWVm8ny+2H2yx9vnFm/nm2D4haJVI6MRduBAREZHI5g8PaUkJTMrvhOWr8fWROhSqrpVb9wdUt2rbwZBN7PYPI9tTVsm+sipyM5LokJYU0cPLJDopXIiIiEjEMsZE3ayEQG/WDYSsF2bJxj08Nr+QBYW7sdY3eX7KkC7cPGUAJ/TpGJI2SHzQDt0i7UhTb0RE4osBThnYOaDayQNy27cxNV4tKObKvyzh4w2+YAHgWpi7uoRL/7SI977YoZ9XEjQKFzHA1HQVu24Au2hLSPn/szZaLk9EJCp53a9uul235RtwxzHk5aZzUgDB4ZqJefWuH2zWWvYcquQnr31GUy9T5bXc/tJKDld5260dEl8ULmJAcnIyAGVlZWFuiTTk/574v0ciIhIdrLVUVnt5dXkxP3x5JXf+YwUvfrKVigBuwr2u5feXjKRnh9Qma66ZlMfXhnVt9/kOzy/ZwhFv8x8+Hqqs5pWC4ppdxEXaRnMuYkBmZiYVFRXs2LEDgPT0dN8YVX1aHhbWWqy1lJWV1X5PMjMzw9wqERE5FvPWlHDHP1Zy4HBV7bFXl3/J/769mt9ePIKzj+/e5HM9jqFzZjL/unkyj80v5JWC4trrjOiVzbWT8rhwdC+ste36s9oYw+LCPQHVLircw9UT89qtLRI/FC5iQE5ODmVlZZSXl1NcXBzu5kgDaWlp5OTkhLsZIiISAK9rWbl1Pzc8X0CV9+hP8ksrq7nlhf/SIS2R8f06Ndnz4HEMHdOS+PnXh/Hjs4awq7SS5ASHThnJtUOhQvEhYKC9Ee04OkvijIZFxQDHcejduzedO3cmJSVFPRYRwBhDSkoKnTt3pnfv3jiO/qmJiEQDj2OYNXd9o8HCz+taHvrP+haHNPlXgkr0OPTokEqnjOTa1wiVE/p0CKyudwdN6pagUM9FjHAch9zcXHJzQ7PyhIiISCzaVVrJB+t3tVi3dNNetuwtp09OWgha1TqutVw5oS9PfLix2Z6J5ASHy8f1Dl3DJKbp41QRERGRGtsPHCbQD/C37T8c0Z/2O8bQq2Mad04f3GzdL74+jOy0JI18kKBQz4WIiIhIjZz0pIBrO0bJDflNUwbQLSuFx+ZvoHDXVytLDu2eyYxpAznruKYnp4scK4ULERGRKFXtuhgMFkuC5nYFRa+OaYzq3YEVW/c3Wze4ayaDu0XPSoAXjO7JRWN68VnxfnYfOkLX7BSGdc8KaO8OkWOhcCEiIhJFXGtxjKGssprZq3awv/wI3bNTOX1YVxI98REwvK7F4xj2lR9hd2klOelJtaswtXWytNe13HBqPjc8X9Bs3fdP7V/7vYgG/nYe1zMb14L/j8kJ4eRyiQ8KFyIiIlHCtRZr4Vdvf8HfFm+pt6tyTnoSN08ZwPUn9QtjC0NjxdZ9PDKvkA/WleBaMAYm5+dy45R8JuW3bWETj2M487hu3H3OUH799upG51/ccfogLhzdq02vEy7GGDzKE9KOFC5ERESihGMMt7ywnDc/3X7Uub1lR/jlv7/gUGU1M6YNDEPr2p+1ljmrdnLT35dTXWc4j7Xw8YbdLCzczf9dMpKLgnDj/92T+3PGsK48u3gzSzbuxWI5sW8OV0/sS//OGW2+vkisUrgQERGJAl7X8t8t+xoNFnU9Mm8DV07oS8e0xKiYbHwsyo94ueMfK+oFi7pcCz959TOmDO4SlPffq2MaPztnWP3X0BwFkWbFx+BMERGRKOdxDM8v3tJi3RGvy4tLW66LNq61/GPZVsqOeJutO+J1+fuS4Lz/xuYjaI6CSPMULkRERKLE+pLSgOrW7TwUc70WjjEs3rgnoNqFhbtj7v2LRAsIFoe4AAAgAElEQVSFCxERkSiRnBDYj+3kxNj88R7ofnUauCQSPrH5v4+IiEiMsdYybWjXgGpPH9oVb4zNDbDWMqp3h4BqTwiwTkSCT+FCREQkSlwxrg+piZ5ma/rkpDF1SJc27/cQiS4f14ekFvbycAxcNbEvbqDdHCISVAoXIiIiUcAYQ4e0RB791ugmh0flZiTx5NUnRtSwoMZWV2rNikvGGDqmJ/GLbwxrtu6uM4fQPTs1aja3E4k1WopWREQkShhjmDK4M+/cejJ/+XgTb6zYxqHKajpnJHPJib24bnI/OqUnRcSKRl7X4nVdXiko5oWlW9m0u4zUJA9nDu/GdZPzWr1XxJUT+tIpI4k/zN3AF9sP1h4f2CWDG6cM4IITegbrLYhIKxirbsOIZYwpGD169OiCgoJwN0VERCKI17W1w56qvS4JNUOFXGsj4hN7r2s5VFnNFX9ezKptB486n+gx3H/JSM4d1bog4LoWxzGs3n6QnQcryM1I5rie2RHz/kWi0ZgxY1i+fPlya+2YtlxHPRciIiJRpu58ioQ6cxAi5cba4xju+MeKRoMFQJXXcufLKxneM5t+ndKPuafFXz+0exaDu2bif9uR8v5F4pnmXIiIiEjQuK5l854y5q0pabauymt5ekFRm4dwOY7RnhZxTCNwIo/ChYiIiASNMfDO5zsC2pNi9qod7d8giTn+ZZaL95WzaOMeVm07AKAVwiKEhkWJiIhI0BhjKD9SHVDt4SPedm6NxKLlW/bxwJx1LKqzY3v/3HS+c3J/rhjfB2uterPCSOFCREREgsbrWvrlBrYSVF5uWju3RmKJtZb5a3fxveeWUeWt30uxcXcZP/3nZxTtKeOnZw8NUwsFNCxKREREgsjjGM4+vhsd0hJbrL1iXJ9W7Xkh0aPh97ctO8d7XcsPX155VLCo64kPN/Jp8f6I2KG+4XyQeBm2pZ4LERERCarkBA8/nD6Yn73+eZM1w3tkcdHoXhGxJ4cEn39Z4E17ynhx6Va27T9MVmoC547sycT8TvWWUw6E17W89dl29pQdabH2uUWb+b9LRral+UGxq7SS55ds5tPiA3gcw0kDcrn0xN6kJXlietiWwoWIiIgE3ZUT+gLw+zlr2V9eVXvcGJg6pAsPXDqq3jK6Ejtca/G6lttfXsEbK7bVO/fC0q2M6JXNX64ZS056UsABw+MYCjbvC6h2WYB17enx+YX835y19XpQ5q4u4ffvruWhb47i9GHdwti69qVwISIiIu3iygl9ufTE3rz12XY27T5EWmICZx3fjb6d0o/5k2uJHo4x/L9XVx4VLPw+LT7At55czNszTgYC/zsQ6Kf94fxb5VrLa8uL+c3sNY2eLzvi5aa//ZdXfzCRYT2yY/LfgD4yEBERkXaTlOBw7sge3HH6YL5/an/65PgmcYf7pspaW29MfLyMh29vrrUU7y3ntf9+2Wzdup2H+Pen2wOeG+F1LeP75QRUO75/YHXtwTGGR+ZtaLbmiNflsfmFYf830F7iLlwYYy42xvzBGPORMeagMcYaY55v5bV6GWP+aozZZoypNMYUGWMeMsZ0DHa7RUREopXHMbiub3lQa6G0oooqrwu0bYJva/gDxcHDVTy9sIj73lrNI/M2sG3f4bC0J9Y4xvCPgq0B7XPycsHWYxoWdcbwbnTNSm6x9pqJeWH5PrrWUrB5L0V7ylusfe+LnZRVBrZkc7SJx2FRPwNGAoeAYmBIay5ijMkHFgJdgDeANcA44FbgTGPMZGvtnmYuISIiEhdca9lZWsGfP9zIKwXFHKyoJtFjOHN4N757Sn9G9OoQknZYa7HAb99ZzVMLiqisdmvP3f/eWs4c3o37Lx1JcoInZj9VDoXtByoCqtsRYJ2fMfDQZSdw3dNLqahyG625/fRBDOmedUzXDRbXWkpKKwOqrXYt+8qPkJ4ce7ficddzAdwODAKygB+04TqP4QsWM6y151trf2ytnQo8CAwG7mtzS0VERKKc17Ws21HK2Q9/xF8XFHGwwvdpbZXX8uan27nwsYW8saL5ITTBYozhvrdW8/gHG+sFCwBrfTuLf/vpZcTwQj4hkZOWFFBdxwDr/BxjmNA/h1dvmMT0YV3rBcDje2bzh8tP4NZpA4/pmsHkGEPnjJZ7VsDXE9PhGN9/tIi9uNQCa+37/l+3dhkwY0x/YDpQBDza4PQ9wPeAq4wxd1pry1rXUhERkehnDHzvuQL21Vkxqq7qmr0LxvXLoWtWCk473dm71rK7tJKnFxY1W7do4x7mr93FqYM6q/eiFbyu5cLRvfjThxtbrL3ghJ64rj2m5YiNMQzpnsUTV5/IvrIjbDtwmOzURHp1TAv7kDbHGE7My6FPThpb9jY/NOprQ7uSEYO9FhCfPRfBMLXmcY61tt5HH9baUmABkAZMCHXDREREIoXXtcxfW9LijVaV1/Lcos3tFizAd+P34idbA7oBfWHpFgWLVvI4hsHdMjltUOdm67plpXBhK/c58X9vOqYnMbxHNr06RsYiAeAbenfTlAHN1iR6DD84NT/sYai9KFy0zuCax3VNnF9f8zgokIsZYwoa+6KV80FEREQigccxzFuzK6Da+WsDq2uLzXsCG0wQaJ00zutaHrliNGPzGl/fpkd2Cn/77niSEmLvNtQYw2Vje/PD6YNpLOukJnp45PLRjOrTISLCUHuIzf6Y9pdd83igifP+46GZoSYiIhKhqr2NT7xt6EiAdW2RkZIYWF2MDlcJFY9jSE3y8PINk/hgbQkvfrKVL/f7hi+dO7IH547qQYLjxOzNNcDNUwdw0ZiePLdoM599eQCPMZw0MJfLxvaO+b9fsf3uwsf/ryWg/i5r7ZhGL+LrvRgdrEaJiIiEkmttwCv3DGvnFX68ruWc47vzTAtzLgDOPr47rrXtOkwr1vmDw0kDO3Pq4C61x+Np88SuWSncdeZXg1Di5e9U7PVHhYa/ZyK7ifNZDepERCSGWGup9rpUe11tvtYMxxguGdOL1ERPi7VXTejbrmPQPY5hXL8cju/Z1I9un+zURC4b27veTaDbSLtidbx8sDUMEvESLICjgkQ8BAtQuGittTWPTc2p8K+D1tScDBERiUL+G8r95VXM/nwHb366ncKSQ/XOSX3pyQn85OzmpxCeN6oHY/vltPuNp9e1PHnNifTPTW/0fFZqAn+9duxRew8U7jrEz/75GZN/M4+xv/oPV/1lCXO+2IHrWn3fRRrQsKjW8S9nO90Y49RdMcoYkwlMBg4Di8PROBERCT6vaymrrOaef63irU+315sjcGLfjtx73nCGdc9q9TLnsezqiXmkJyXw4H/WUVyzEzZAVkoCV03M487TBx3zkqSt4XEMuRnJvDXjZF5bXsxLy7aybf9hslOTOHdUD64c34eO6Un1PmH++5LN/Oz1z6mbIXatr+Sj9bs5bXBn/nz1iRiXdm+7SLRQuGiGMSYRyAeqrLWF/uPW2kJjzBx8e13cBPyhztPuBdKBP2mPCxGR2GCtpbLKy2VPLGL19tKjzi/bvI9LHl/EP74/kaHds+Jq6EegLhzdkwtG92Thht0U10zunTq4C8mJHqy1IQtl/snGl4/rw7cm9K09bmuGt/nb4XUtK7bu4+7XP6epkW/z1+7iV2+t5t5zh7d7u0WiRdyFC2PM+cD5Nb/tVvM40RjzdM2vd1trf1jz657AamAzkNfgUjcCC4FZxphpNXXjgSn4hkPd3R7tFxGR8PjLgk2NBgu/8iNeZv5rFa/8YFIIWxU9jDEYYGJ+bu2NfILHqT0Xag17Ghq2weMYnvxoU5PBwu8fn2zlrjMGHzWUSiRexeOci1HANTVfZ9Qc61/n2MWBXKSmJ+NE4Gl8oeJOfL0cs4CJ1to9QW21iIiEjQVeWLKlxbplm/exfmepJnk3w+MYEjxObbCIVNVel/e+2Nli3eEqL/9ZvZNqt/2X0hWJBnEXs621M4GZAdYW8dWyso2d3wpcF4x2iYhI5CqtqGbbgYqAaj/78gD9ctNxPBoaFc2OeF2qA5ysXVZZ3c6tEYkekf2xgYiISARIOIY5FIkR/om8BCYtKYHOGckB1eblpmOa/ixSJK7of0AREZEWpCcncELvDi3WJXkcThqQqwndMeLSsb1brOnbKY1J+fqei/gpXIiIiLTA61qumZTXYt3Zx3enY3qSlqONAa61fPukfnTPTmmyxhi464zBtRPURUThQkREpEUex3D+CT35ZjOfZA/tnskvzxuuTdVihGMM2amJvPz9iYzodfSu3h3TErn/kpGcM6KHwqRIHXE3oVtERKQ1XGv5zUUjmNC/E88sLOK/W/cD0CM7hcvH9+Hbk/uRnOjR8JgY4nEM3Tuk8q+bT2LF1v3MXb2TI9UuQ7tncc6I7ppfI9IIhQsREZEA+Hdt/sbIHpx/Qk8OVVZT7XXJSk2sncqrT7Bjjz8sHt8zm1E18268rlWIFGmCwoWIiMgx8N9UZmjTtLhSN0woWIg0Tf15IiIiIiISFAoXIiIiIiISFAoXIiIiIiISFAoXIiIiIiISFAoXIiIiIiISFAoXIiIiIiISFAoXIiIiIiISFAoXIiIiIiISFAoXIiIiIiISFAoXIiIiIiISFAoXIiIiccy1lmqvS7XrYq0Nd3Nihtet/2dprdWfr8SFhHA3QERERELP61o8jqF472FWbN2HYwwT8juRm5Fce06Onde1GGDemp38fclWNu0+RHpyAmce141vje9Lx7REjNGfrcQuhQsREZE441rLlr3l/OKNz/l4w278H6gnOIYzj+vGL887juzURAWMY+R1LUeqXa5/+hMWbdxT79yqbQd5fH4hT1x9IpPyOylgSMzSsCgREZE44nUtW/aUc+FjC/ho/VfBAqDatfz70+1c9MeFHKqoxtUwnmPicQx3vbLyqGDhV3bEy3efXcaOgxW4rv5sJTYpXIiIiMQRj2P41VtfsK+8qsmaTbvLeOT9DTj6dD1grrVsP3CYtz7b3mxd+REvzywswlGvkMQohQsREZE4svNgBfPWlLRY93LBVqq8bghaFBsM8ObK7QTSIfHmyuYDiEg0U7gQERGJE65rWbXtYEA3wPvLq9i+/3D7NypGGGM4WNF0b1BdBw8HVicSjRQuRERE4oXxTdoOVIJHtwmBcl1Lzw6pAdX2ygmsTiQa6X8NERGROOEYw4l5HclIbnmxyPzO6fQI8GZZwHEM547sQXqSp8XaS0/srcnyErMULkREROJIWlICF43p1WLdlRP66gb4GKUnJ/D9U/ObremXm85lY3trsrzELIULERGROOK6lh+fOYQTendosuas47pxzcQ83QC3woxpA5kxbQDJCUffYo3slc2L35tAckLLvRsi0Uqb6ImIiMQRxzEkJTi8+L0JPLWgiL8t3czWvb6J20O7Z3LVhDy+Oa43qNOi1e44fTDXT+7HP5YVs2n3IdKSEjjn+O6M7ttRu59LzFO4EBERiTMex+AYh++f2p8bTstnX/kRPMaQlZqI61pfj4Xuf9skOzWR757cr3Ynbm/NEl0KFhLrFC5ERETikKkz5KljWlLtr7W5W3CYBkPKFCokXmjOhYiIiIhIGLmNbD7jDWRDmgikngsRERERkTDauLuMZxYWsWTTHqyF0X07cs3EvgzrkY219qiesEimcCEiIiIiEiaPf1DIb2evoe7Kz+tLDvHSJ1v5/in9+cnZQ8PXuFaI22FRxphexpi/GmO2GWMqjTFFxpiHjDEdj+Ea840xtpmvlPZ8DyIiIiISnbyuZfbnO/jNO/WDRV1/+nAjf1+yOar2nInLngtjTD6wEOgCvAGsAcYBtwJnGmMmW2v3HMMl723ieHWbGioiIiIiMcnjGJ74sLDFuj9/tIkrxvcNQYuCIy7DBfAYvmAxw1r7B/9BY8wDwO3AfcANgV7MWjsz2A0UERERkdi1ff9hlm/Z32Ldpt1lfPblAY7vmR2CVrVd3A2LMsb0B6YDRcCjDU7fA5QBVxlj0kPcNBERERGJEwcqqgKu3V9+JGqGRsVjz8XUmsc51lq37glrbakxZgG+8DEBmBvIBY0xlwH9gCPAamCetbYyeE2WmON6wfHUP2ZdMHGX90VEROJSt6wUPI4JaMnZPjlpvs0to0A8hovBNY/rmji/Hl+4GESA4QJ4scHvS4wxN1lrXwnkycaYgiZODQnw9SVaWAvGwK618Mmf4csCMB7odwqM/Q506B3uFoqIiEgIdEhL4mtDu/Duqp3N1o3N60jfTtEzoCYePyb1D1g70MR5//EOAVzrDeAbQC8gFV8Y+N+a575kjDmrDe2UWGQMvPcL+ONEWPZX2L4Sti2HBQ/BrJGw7Klwt1BEJCANh2hE64ZfIuHidS0zpg0kOaHp23GPY7j99EFR9e8rHnsuWuLvc2rxu2itfbDBobXAT40x24A/AL8G3gngOmMabYivR2N0S8+XKGEtfPIkLHi48fOuF966HTr2hX6nHj1sSkQkAvg39Cree5jnl2xm/c5DpCY5fG1oV74xsgceY3Cc6Bi+IRJOHscwrHsWT183lhkvrmBXaf0R9R3SEvnthSOYlJ8bpha2TjyGC3/PRFNT7rMa1LXGk8CDwChjTKa1trQN15KYYWHhH1oosbBgFuRPbb5OROKGrekhqLtDr2ttWMZf+4PFzH+t4plFRfXW5n/7M996/X+9dizDemRFzfhwkXAyxjC+fycW/WQq73y2g8Ub/Tt0d+DckT1J9ETfv6N4DBdrax4HNXF+YM1jU3MyWmStrTDGlAIdgXRA4SLeWQtbP4H9m1uu3TQfyvdCWk67N0tEIt+esiP8bfFmFm3cg2thVO8OXD2hL71y0kLeFmMMD7y3lqcXFjV6vqS0kiv/soR3bzuFzpnJChgiAXCMwTGGs4/vzjdG9gB8Q6Y8UdoDGI/h4v2ax+nGGKfuilHGmExgMnAYWNzaFzDGDMYXLEqB3W1oq8QK68LhAPdltBYO71O4EBFeLSjmJ699xhHvV4sbLt20lz9/tJE7Th/ELVMHNvPs4Cs/Us1fPy5qtmZ/eRXPLCzirjO1JonIsagbJqI1WEAcTui21hYCc4A84KYGp+/F19PwrLW2zH/QGDPEGFPvf0ljTH9jTM+G1zfG5AL+WbkvWmu1S7f4lpjNDnAlKE8iZHRp3/aISETzupZFhbv5f6+srBcs/KyF++es4+VlW0O29r3XdXnn8x0cqmz5x9orBcUhaJGIRKJ47LkAuBFYCMwyxkzDtzfFeGAKvuFQdzeoX13zWDdGngI8aYz5ACgE9gJ9gLPxzedYBtzVXm9Aoowx0O146DYCdnzafO2Qb0ByZmjaJSIRyeMY/ji/kJYWiPnj/EIuOTE0S1gbY9h1MLAtnHYd0lZPIvEq7nouoLb34kTgaXyh4k4gH5gFTLTWBjJ+pQB4HugCXFRzjTOBz4AZwGRrbct7ukv8cF047Se+oNGUxFQ45U7fylEiErf2lR3how0tj6rduLuMlcX7ayd9tydroVNGUkC1uenJ7dwaEYlU8dpzgbV2K3BdgLVH3Q1aaz8Drg1ysySWOQ4MORvOexTevguOHKp/Pj0XLn4KugxvPoCISMzbV36EQPPC7tJKXAvtvaiMx/FNOL3nX6soP9L8ByAXjD5q1LCIxIm4DRciYTPqWzDsfFjxd/hymW+H7v6nwvALfXtbKFiIxL3cjGQ8jglo46xu2SmEau5nenICV0/sy+MfbGyyJislgesm54VtuVwRCS+FC5FwSEqHsd+Gcd/1/d71atM8EamVlZrI1CFdeO+Lnc3WDe2eyfAeTW3bFHzWWu46Ywj7yqt46ZOtR53vmJbIk9eMpWtWioKFSJxSuBAJF1NnypOChYjU4XUtN08ZwAdrdzW6WpTfjGkDcV0bsh2xjTFYLL+9aATfntyP5xZvZkPJIZITHaYP68aFo3uSlOAoWIjEMYULERGRCONxDCN7d+CxK0dz+4srKG2w/GtygsPPvz6Ms47rHvK2+YNDfpcM/uf842qPhzLkiEjkUrgQERGJUNOGdGHp3V/j1eXFLCrcg9daRvXuwOVj+5CdlhjWtjXc5EvBQkRA4UJERCRiGWNITfJwxbg+XDmhL0BIlp2V0Kp2XQwGsHicuNwlQGKIwoWIiEiEq9srYDSfISZYazHGUH6kmve+2Mn+8iq6Z6cwdWgXEhQwJIopXIiIiIiEkFuzxPDvZq/hucWbOVRnTk3njGRumTaAqyfmhal1Im2jcCEiIiISQo5juP2lFfzzv18edW7XoUp+8cYqDlVUc+OUAWFoXeP8E/b3lx/B61o6pifh70NTb5rUpXAhIiIiEiJe17KyeH+jwaKuh/6znivG9yE7NTGsN++ua8HAS8u28tyizXyx/SAA3bNTuHxcH75zcj+SEzxHTfCX+KVBfSIiIiIh4nEMf1u8ucW6I16XF5cevVFhyBm47cUV/OS1z2qDBcD2AxU88N46Lnl8ERVV3oB2k5f4oHAhIiIiEkJrdpQGVLd2Z2nYey1eXLqVf63c1mTNqm0H+eW/v1DPhdRSuBAREREJoeSEwG6/kgKsay+OY3h2UVGLdW+s+JLSiqp2b49EB4ULERERkRCx1jJ1SNeAaqcN6RLW4UZ7DlUG1MtSUeWyeOMeDY0SQOFCREREJKS+Nb4PKYnN34L16pjK14Z2Detwo+pjCAvVrsWicCFBChfGmBOMMbcbY24xxgxupu48Y8xfg/GaIiIiItHGGEOHtEQeuXx0k8OjctKT+PPVJ4b9Vj03I5nOmckt1hkDx/fIxqMlaYUghAtjzO+BZcDvgYeAVcaYWcaYxEbKRwHXtPU1RURERKKVMYavDevKm7ecxKUn9iY10QNAx7REvndKf9659WQGd80M+yRpj2P45tjeLdadMrAzvXLStN+FAG3c58IYcwFwB3AQeAmoAi4CbgKON8acY60tb3MrRURERGJMfucMfnfxCH538Qgqq70kJ/hChmstTgTcqLvW8v1T8nl31Q7W7TzUaE12aiI///owvK4NexiSyNDWnosfABXAeGvt9621NwODgFeBU4F/G2NS2vgaIiIiIlHD61pca2sfm1L3ZtwfLICICBbga0dqkoeXb5jEhaN71hvGZQycPDCXV38wif6d0xUspFZbd+geDbxmrV3rP2CtLQUuNcY8ANwGvGmM+bq1trKNryUiIiISsay1WAvz1uyksKSM5ESH6cO70bNDatR+su9xDBnJCTxw6Sh+8fVhLCzcg+taRvTKpk+ndLxuZPSySORoa7jIABrdZtJae4cxxgvcCbxujDmvja8lIiIiErHe+mw79721mu0HKmqP/c+/v+D0YV353cUjyUhOiNqAAb4hUGcM74blq0ARje9H2ldbh0VtB7o1ddJa+/+Ah4Ez8A2VSmrj64mIiIhEFGst//50G7e88N96wQLAtfDuqp1c/sRiqrxumFoYHMYYPI4hwXHUWyFNamu4WI1vbkWTrLW3A38EzgFubePriYiIiEQUa+G+t1bTzPQKvth+kJc+2drsHAyRWNDWcPEO0N8Yc3JzRdbam4C/AGltfD0RERGRiOF1LfPWlBzVY9GYvy/Zok/8Jea1dc7Fq0APoFNLhdba7xpjioG8Nr6miIiISEQwBjaUNL5Ma0MbdgVWJxLN2hQurLXbgJ8cQ/29xpi2BhoRERGRiGAtpCR5Wi4EUprYkVsklgTtb7kx5omW9rQwxuQBHwfrNUVERETCyQDTh3UlkEWTpg9vcg0ckZgRzAj9HWCpMWZIYyeNMRcD/wXGBvE1RURERMLGcQw9OqRy+rDmg4Nj4LrJeXhdTeiW2BbMcHEfMAxYZoy5zn/QGJNkjHkMeAnwAhcE8TVFREREwsrrWn538QiGdc9q9LwxMPPc4Yzo1UH7QkjMC9r8B2vtz40x84HngSeNMVPx7XHxJDACWABcbq0tDtZrioiIiISbfxfr126cxEufbOXvS7awYdchUhM9fG1oV66bnMfI3h3C3UyRkAjq5Gpr7VxjzEjgOeCKmi8X+BUw01ob3bvHiIiIiDTC4xg8joerJvblmkl59c5pKJTEk/ZYuekQsAvfHCeAA8CHChYiIiIS6xrbx0JDoSSeBHVNtJpei+XA5cC7wA1AEjDbGHOfMUZrsImIiIiIxKhgLkV7E7AI6A/81Fp7lrX2CWAM8CnwY+AjY0yfYL2miIiIiIhEjmD2JPwBKAFOtdb+1n/QWrsemAA8BkzEtxxtWBljehlj/mqM2WaMqTTGFBljHjLGdDzG6+TUPK+o5jrbaq7bq73aLiIiIiISqYIZLt4ATrDWLmp4wlp7xFp7C3BhEF+vVYwx+UABcB2wFHgQ2AjcCiwyxnQK8Dqd8PXU3AoU1lxnac11C4wx/YPfehERERGRyBW0cGGtvcBau6+FmteBUcF6zVZ6DOgCzLDWnm+t/bG1diq+cDAY334dgfg1MAh40Fo7reY65+MLG11qXkdEREREJG6EfIK1tXZrqF/Tr6Y3YTpQBDza4PQ9QBlwlTEmvYXrpANX1dTf0+D0IzXXP0O9FyIiIiIST+Jt9aapNY9zGi6Na60txbfRXxq+OSLNmQikAgtqnlf3Oi4wp+a3U9rcYhERERGRKNEe+1xEssE1j+uaOL8eX8/GIGBuG69DzXVaZIwpaOLUkECeLyISE1wvOB6oOgzeI5CcCcYB6/oeRUQk4sVbuMiueTzQxHn/8Q4huo6IiPht+A8s/iNsmg/WQkYXGH01TLwZkrN8wUNERCJavIWLlvi30LShvI61dkyjF/H1aIxuY1tERCLf/N/A/P+tf+xQCXz4e/jsFbj+XUjvrIAhIhLh4q2f2d+jkN3E+awGde19HRGR+OZ6oWjB0cGirn1F8M8bFCxERKJAvIWLtTWPTc2FGFjz2NRcimBfR0QkvjkeWPJ4y3Ub34fdG3xhREREIla8hYv3ax6nG1N/dgvSHVoAACAASURBVKAxJhOYDBwGFrdwncU1dZNrnlf3Og6+SeF1X09ERJqy4T+B1a2f3b7tEBGRNourcGGtLcS3TGwecFOD0/cC6cCz1toy/0FjzBBjTL1Vm6y1h4DnaupnNrjOzTXXf9dauzGIzRcRiU1uVWB13mowpuU6EREJm3ic0H0jsBCYZYyZBqwGxuPbk2IdcHeD+tU1jw1/ov2U/9/efYfZVZWLH/++Z9IDSSiB0EMPvUWqQgBRQCkKKCgiIIjlXvvvWrFcr9fGtWG71wihCooCFkBqQERQEor0GloiLRBKSJuzfn+sc5jJZMqZzJk57ft5nvPsOXuvs2btlZ05+92rwTTgUxGxI/B3YCvgMOAZVgxeJEndWXs7mDu773STtndaWkmqcy33F7rUejEVmEEOKj4NbAr8CNgjpfR8hfk8T15M70fAZqV8dgPOBHYp/R5JUm+K7TD1xL7Trb4JbLYfFFrxmZgkNY6W/CudUnoCOKHCtD22waeU5gMfL70kSf1VaIMdjoY7L4A5N/aQZhgcfJqtFpLUAPwrLUmqrSjAsb+FqR+A4WOWPzZpezj2d7DZ/gYWktQAWrLlQpJUR8rrV7z9e3DA1/LsUUsXwsStYb2dcouFJKkhGFxIkmqvHGCMXBW2OqS0s9Qr1RYLSWoYBheSpPrioG1Jalg+DpIkSZJUFQYXkiRJkqrC4EKSJEkVaS+m138uptRLSrUqO7ZKkiSpV+3FRFshuPmR57n2vmdYvKzIlEmr8o6d1mPsSG8n1cGrQZIkST1qLyYee/5VPnTuLB54+pXljn3zsnv57EFTOG6PybUpnOqOwYUkSZK6VUyJFxYu4d3/dzPPvrx4heOvLmnny5fezZgRbRyx8/pERA1KqXrimAtJkiR1qxDBGTc+2m1g0dn3rnwAh2AIDC4kSZLUiwv/8USfaeYuWMTMB56lvVgcghKpnhlcSM2u2J63C5+HZ+6BV55efr8kST1YvLSd519dUlHaJ19YCNgtqtU55kJqdk/Nghu+Cw9dDan0RGnym+CNn4TN9q9t2SRJdW34sAIjhxVYvKzvFolxo4aTSBhgtDZbLqRmlRLc9yeYcTA8eGVHYAEw5y9w3hEw++zalU+SVPcKERy47aQ+040Z0cYBW6/NsIK3lq3OK0BqVktfg4tPgfal3R9PCf74ydxNylF4kqRutBcTJ79pE4YVem+NOGbXDV3vQoDBhdScUhHuOB8Wv9x7uuIy+McvwakDJUndaCsE2643ntOO2oHhbd1/V7xl67X5/EFTXLFbgGMupOYUBZhzY2Vp5/xlcMsiSWp4h++0HlMnr8Y5f3uMa+97hiXtRbZce1WO3X0j9t5iIsViouCDKmFwITWvVOF0gD5pkiRVYN0Jo/n8wVvx+YO3en1fezF/hxT66Dal1mG3KKkZpSKs/4bK0m6w6+CWRZLUFLprmWgzqFAXBhdSUwrY+TgYNqqPZAFvOKnyVg5JkqReGFxIzSgCRo2Hg0/rfbD2fl+GCRvmMRqSJEkD5JgLqZnt/D4YuyZc/22Ye1vH/olT8iJ6Oxxdu7JJkqSmY3AhNbvN3wJbHgTP3Acvz4Uxa8I629sVqlUV26HQ1vG+PKDfWV4kSVVgcCE1u/KN5FpTYOIWHV2g7ArVWortOYB4+DqYdSY89wCMWAW2ejvscgKMnuA1IUkaMIMLqZV489iaiu1QXAq/eg88fM3yx+bOhr98D44+DzbexxYMSdKAeKchSc2u0AZ//OSKgUXZklfggvfAgiftLidJGhCDC0lqZinBK0/Dnb/uPd2SV+Hv/2frliRpQPwWkaSmluCu30FxWd9J77po8IsjSWpqBheS1MyiAIterCztaxWmkySpBwYXktTMiu0wfv3K0laaTpKkHhhcSFIzK7TBNkfAyFX7TrvzcQ7oliQNiMGFJDW7EWPyiuy9mbAR7HI84FS0kqSVZ3AhSa3gTZ+Gff4D2kaseGzSdnD8H2H4GNe5kCQNSEsuohcRewJfAnYHRgEPAWcAp6eU2vuRT+rl8C0ppd0HVFBJqqZ9vwi7fghuOweefyi3aGx1KEx+Yx6bUV7NXZKkldRywUVEHAb8FlgEXAjMBw4Bvg/sBRzVzywfA2Z0s//JlS+lGkpKUI5Jo80nv6pvY1aHvT7ecZ0WS9eugYUkqQpaKriIiHHAL4B2YFpK6dbS/lOBa4EjI+LolNIF/ch2Tkrpq1UvrOpf+UnvK0/DE3/PN2sb7gFj1/QpsOpX1+DX61SSVEUtFVwARwITgbPLgQVASmlRRHwJuAb4MNCf4EKtqNgOL8+Dyz8LD1ze8fS3bThMOQQO/g6MXr37G7dUzGsPLHkViNw1xWBEkiQ1gVYLLvYrba/o5tgNwEJgz4gYmVJaXGGeEyLiRGASsACYlVK6eeBFVd0qtufWiulvzgFGZ+1L4e7fwdzZcNK1MHpCR9BQbAcSzDoL/jEdnrkn719nB3jDybDje/NxgwxJktSgWi242LK0faDrgZTSsoh4FNgG2AS4t8I8dwB+2XlHRNwBvC+l9M9KMoiIWT0cmlJhGTSUCm1w9VdWDCw6e2EOzPxveNv/5PepmF/nvxsevmb5tPPugN//Gzz4Zzjq7DyGw3EbkiSpAbXaVLTjS9sFPRwv759QYX7fIw8CnwisCrwBuIgccFwbEeutZDlVz157Ee65tO90d14IS1/LP0cBrv/2ioFFZ/f+Af724+qUUZIkqQYaLriIiDkRkfrxOrc/2Ze2vU0x+7qU0qdTSjellJ5LKb2SUro1pXQUeTaqNYHPVJjPLt29gPv6UXYNhZTg6btgWQW95ha/DM8/nH9uXwqzzuz7M/+YToWXX2sptvf+XpIk1YVG7Bb1MHka2UrN7fRzuWVifHcJgXFd0q2snwNHAHsPMB/Vo0I//tuU0869DV59ru/0Lz4Gzz0IE7fsO213ug4Mb/SB4uXB7889ALeeAfMfgZHjYOvDYKtDcppGPj9JkppMwwUXKaX9B/Dx+4GpwBbAcuMcImIYsDGwDHhkAL8D4NnSduwA81G9iYB1d4bRq8FrL/Sedtx6sObmubWjfWnlv6N9Sf/LVWzPa23c9Vu4+2JYtADGbwA7vx82flNjBhnlcSoXnwJ3/nr5Y3f/DlabDO+9CFbfpPHOTZKkJtVw3aIG6NrS9sBuju0NjAFu6sdMUT0pr8w90CBF9WjYCNj5uL7TTT0hP3WPgLW2ytPU9mXE2Hyz3B/F9vxE/0c7w8Ufggf+DI/fDP/8DZz1djj7cGhf3HhdiaIAl/2/FQOLshfmwNmHlqb0lSRJ9aDVgouLgOeAoyNianlnRIwC/qv09medPxARYyJiSkRs2GX/zhGxQstERGwPfKP0tj/jPdQoikXY9wuw4e49p9l0f9jrEx3vx6ye17/oy3ZH5QCjP5YtyjfZC57o/vgj18HvPthYT/dTglefhdln957upbkw+6zcwiFJkmqupYKLlNJLwMlAGzAzIqZHxHeA24E9yMHHhV0+tit5WtqudzkfA+ZFxCURcXpEnBYRfwRmA2uQVwL/1eCdjWqmUIDCcDjuD/Dmr+XuOWVrbAYHfhPe++uOVgvIrQZv/gqMWaPnfMetC9M+378WhlSE28/LN9m9ufcPeXB5sYFuwu+4AIrLKkj3q1zXkiSp5hpuzMVApZQuiYh9gC+SB12PAh4CPgX8KKVU6VQ9l5AHgG9PXpxvFPA8cDnwi5TS76tddtWRQlu+od3rY/n16nNAwCoTOwYhd00/YUP4wFXwp0/BIzM7jkXklo5DfghjJ/avhSEKecrbStxxPux3auV511JEXqiwEpWmkyRJg67lgguAlNJfgYMrTDuTjilqO++/hBxgqFVF8Pqlscpanfb38BQ9CrmV47hLcyvCEzfnz09+Yw48VnbQ9cL5laV79fn+510rKcHYtfpOBzkgkyRJdaElgwupZsrBw+obw2ob5SUtCoXlj/XXKmvnAd19WXXtlcu/VrZ/N1zztb67Ru1wTPetRZIkacj5bSzVQhTyGhhtwwZ2U1xshx3fU9nv2/HYxhlzEZGDob7ObZW183S7KzYuSpKkGjC4kBpZoQ22fxessWnv6bZ/N0zYoKOVpBGkIrzte7D14d0fH7ceHPd7GLlqx8B5SZJUU3aLkhpdYTi8/w9w3lHw9N0rHt/2CDj09MbrOhSF3CDxrrNg7u0w60yY/2ieqnebw2Gbd5ZagBpoil1JkpqcwYXU6AptsMok+PBN8ODVcNdFeYXuCRvCLsfnBfyKDRZYlJVbWiZtl2fTKmvEFcclSWoBBhdSMyjfaG+6L2z+5o795TUzGqk7VHe6BhIGFpIk1aUGv+OQtBxvwiVJUg0ZXEiSJEmqCoMLSZIkSVVhcCFJkiSpKgwuJEmSJFWFs0VJklZUXhelfSk8fQ8Ul8LEKTByFacCliT1yOBCkrS8Ynt+3fBtmHUWvPps3j98DGx3FOx/Koxe3QBDkrQCgwupWlKCiHxTtuDJ/PP49cnLTJPfS/Uupdxqce47YM6Nyx9buhBmnwWPXAcnXQ1j1jTAkCQtx+BCqoZiOyxbDH/9fn7S+8rTef+EDWHqibDHv+UuJo14I1Ys5kX4li2GZ+8DInePGTai45iaSIKbf7piYNHZi4/D5Z+Do84cumJJkhqCwYU0UMV2WLYIzjwY5t2+/LEXH4ervwoPXwfH/hZSobFaMIrt+Wn1zG/CbefBohfz/tGrwU7vg2mfh2EjGzNoUs9uPaPvNPf9IXeXGrNmY13TkqRB5SNHaaAKbXDN11YMLDp79Hr46w+GrkzVUA4szngr/O0nHYEFwGsvwE0/gjMPyoFVsVi7cqq6Xn4aXpjTd7r2pfDELZDaB71IkqTGYXAhDdTShXD7r/pOd+uZuS97oyi0wbVfh6fv7jnNvNth5rfsGiVJkgCDC2lgUoKnZsPil/pO+9JT8MKjg1+maln6WmVB0+3nwrIlg18eZSnlV1mxyi0Hq64Nq03uO13bcNhgNwi7xEmSOhhcSAPVn5u7at8IDqZ//bOyoGnhfHju/sEvj7LXXoCbTofffwyu+Dw8+Y+8v5rX1tQT+04z5e0wdqLjLSRJy3FAtzQQEbDO9nlQ87LFvacdvRqstvHQlKsqUt9JXk/aj7Raedd9A278AbR3aim6+aew7k5w9HmwyqQqDK4P2P0j8OCVPc8YNWFDOOjbLqYnSVqBLRfSQI1eDbY+rO90O743T9/aKNbeDkau2ne60avBxC0Hvzyt7obvwvXfWT6wKJt7G8x4ex5cP1ARedrkYy+GvT+TWyfKho+BXY6Hk691jQtJUrdsuZAGqtgOB3wdHrspL57XnYlbwj7/kQd0xxDG9J27ypRvGis1Ygxs/274x/Te0+34ntxyo8Gz5NXcYtGb+Y/A7LNht1MGfo0V2vL1st+psM/n8qD+4rK8vsnIVWyxkCT1yJYLaaAKbfnp7klXw7ZH5IGuZcNH5/UgTrwSRqwydIFFeVaqx26CW36WA4QXHsv7Ku2bX2zPN5e9tUqsvU1e68KpaAdPsR3u+i0seaXvtLPPrt41Vs6nbTisuyOsPzUHFmBgIUnqkS0XUjUU2mDsWnDkGfDqczB3NhCwwa4wavzQP+mdcyP86VPw3IMd+674LGy6Hxz+cxizRt/lKbTlblEnXpmnpL3zQlj8cj42ajzscHQOPoaPGZqpaAfSCtPICm15McZKLHhicMsiSVIfDC6kainfrI9dEzbdH4iOm+6hCiyK7fD4zXDuESv2zU8JHromL4p38nUwalzfN+jlAONt/5O7fj19V94/abvcKjMUQVO5K9ljN8G/7oDCcNj8LbD6xq3RPScVczBXiUrTSZI0SAwupMFQqxveQhv8+QvdD/otm/8I3PwT2PeLlecJeQzGBrt2f2wwVaMVptFtewRc9eW+F2Hc9sihKY8kST1okX4FUgsoFmHuHXnV7L7MOqv+VwsvtsOcv+ZWmM6BBSzfCrP45fo/l4GIAoxbF7Z9Z+/pRoyFXU92/IskqaYMLqRqKC4rvWp5Y1eEZ+6uLOkrT8OiBYNbnIHqTytMs4+/KLbDoafDxnt3f3zEKnD0+TBuvaEZ/yJJUg/sFiWtrGIx38gtWQgPX5OnC11zC1hv545jQ2346MrT1vP0scViXiG80laYaZ9v7gCj0AaMhPddCg/+GWadCc8/nFsrtjoEdjkRxqzmatmSpJozuJBWRrE9d8X585fy9J/lWZQA1t4WDvhP2Gz/oS1TYVgeSD58NCx9rfe0G++TZ3mqW0V49p7KkpZbYUavNrhFqrXyuJLN3wJbHtSxv7w6uoGFJKkONPGjPmkwBVxwDPztJ8sHFpBnVDr/KLjvso4bv6Eyahxsf3Tf6Xb7UOXrXdTKsFH9SFvHrTDV1nXweoSBhSSpbhhcSP1VbId7L4UHr+o9zWWfGfqBxsV2OPCbsOEePad502dgysH1PcNS51aYvtR9K4wkSa3D4ELqr0Ib3HpG3+leegoeuGJoWwgKbdA2At7/BzjkhzBp+zwWoW04bHkwHHcp7H/q0LeorIxmaoWRJKlFtNSYi4gYDnwE2BHYCdgaGA6cnFKavpJ57gl8CdgdGAU8BJwBnJ5S8o6nWc27o/J0U942uGXpqtAGtMFOx8Eux3csQgcdN+GN0I2m3Arz7H3w+N+6T1NuhZEkSXWhpYILYCzwg9LPTwP/AjZY2cwi4jDgt8Ai4EJgPnAI8H1gL+CogRRWdaxQ4X+dtuG5laAWN/Pl2ao6z6JUz12huiq0AaVWmNvPg3/8Mo9nKbTBZgfAbqfAJtNqV7+SJGkFrRZcLAQOBm5PKc2LiK8CX1mZjCJiHPALoB2YllK6tbT/VOBa4MiIODqldEFVSq76stmb4c4L+063+VtLLQcNdFNfT5qlFUaSpBbRUmMuUkpLUkqXp5TmVSG7I4GJwAXlwKL0OxaRu0kBfLgKv0f1ptie+/n3Zf2psM72jdVaUK8avRVGkqQW0VLBRZXtV9pe0c2xG8itJHtGRAvNkdkiCm15obw3f63nNOPWhSPOcKCxJElqKQYXK2/L0vaBrgdSSsuAR8ndzjYZykJpCL3xE/De3+YpU8tdc0avBnt+DD54A4xf36frkiSppbTamItqGl/aLujheHn/hL4yiohZPRya0t9CaYhtui9s/mZYshCWvQajJuSAovPYAEmSpBbRcHc/ETEnIlI/XufWqqilbQMsKKCVVm6ZGDEGxqzR8d7AQpIktaBGbLl4mDz1a6XmDlI5yi0T43s4Pq5Luh6llHbpbn+pRWPn/hdNkiRJGnoNF1yklPavdRlK7gemAlsAy3VriohhwMbAMuCRoS+aJEmSNPTsu7Hyri1tD+zm2N7AGOCmlNLioSuSJEmSVDsGF32IiPERMSUi1uly6CLgOeDoiJjaKf0o4L9Kb382RMWUJEmSaq7hukUNVER8jo5ZmHYsbU+IiDeWfr4xpTS900feAZwJnAUcX96ZUnopIk4mBxkzI+ICYD5wKHma2ouACpZwliRJkppDywUX5G5M+3TZt2fpVTadCqSULomIfYAvAkcAo4CHgE8BP0opOVOUJEmSWkbLBRcppWn9TD8DmNHL8b8CBw+oUJIkSVITcMyFJEmSpKowuJAkSZJUFQYXkiRJkqrC4EKSJElSVRhcSJIkSaqKlpstSpLUQIrtpR8SRCG/JEl1y+BCklR/UjEHEk/+A/51JxSGw+YHwPj1c8BRaKt1CSVJ3TC4kCTVn8f+Bpd9Gp65t2NfFGCLt8KhP4HREwwwJKkO2b4sSaofxfYcWJxz+PKBBeTWjPsvhzPeAktehZRqU0ZJUo8MLiRJ9aPQBld8DtqX9Jzm+Yfg5p8OXZkkSRUzuJAk1YdiEebdAfNu7zvtrBmALReSVG8MLiRJdaIIT99VWdKX58HiVwa3OJKkfjO4kCTVj2GjKk/b5pwkklRvDC4kSfUh2mDT/WDYyL7TTn4TDB8z+GWSJPWLwYUkqT5EwOjVYLuj+k6724c6LbAnSaoXBheSpPpRbIeDvgPr7dJzmr0+AVu93XUuJKkO2WFVklQ/Cm153MUJl+cZoW49A569Ly+gt+l+ucVi8wPyGhcRtS6tJKkLgwtJUn0ptOXXrifDbqdAcVkOLqLQ0RXKwEKS6pLBhaTGUr7RTAkIKNi7s2lF6d+20Omryq5QklTXDC60olTsdPOWOr7gpVpKKV+bD1wJz90Pw0fDlm+DCRvkp9nedEqSVHMGF8rK/ZcXvwz//A0seBJGTYDtjoBx63nzptq7+2K48kvw0lMd+674PGx5EBz2Uxi5qteoJEk1ZnChjsDiuv+Gm06HpQs7jl39Fdj2nXDoj4ER3rxp6KUEd/8OfvuBUmta52NFuO9P8MLb4KSroTC6NmWUJEmAU9EKcmBx9dfg+m8vH1hAvnn750Xwq6MdQKnaSMXcYtE1sOjs6btg9tk5rSRJqhmDi1aXErz6LPzt9N7TPTITHrzaRas0tIrL4MEr4aW5faeddWbjjg8qdgmK/H8mSWpQdotSfuLbvrTvdLPOhC3eMvjlkcqiAM89UFnaStPVk3KXxPkPw62/hGfLA9UPgu3eBW3DGzdgkiS1JIOLVhcBzz9cWdr5FaaTqiWlvKBaJYY12HiL8mxsf/x0Diw6u+9PuaviMRfCejsZYEiSGobfWIIRY6ubTqqWCJjytspurqccPPjlqaYIuPYbKwYWZa8+C+e+E17+14rdpiRJqlMGF62u2A5bH15Z2q0Oc8CshlYUYPz6OcDoK91uH2mssQpLXoVbft57mkUvwi3/60KBkqSG4TdWqyu0weS9YJ0de083chzs8n7AGaM0xIrteSrkSdt3fzwKcPB3c/ehRpkqubgM7rkElrzSd9o7fjX45ZEkqUoMLpRv3o45H1bfpPvjI1eFY34Fo8Y7Ha2GXqEtX4MnXQUHnwZrbwOFYXnf9u+Gk66FN5xU61L2TxTglWcqS/vqs4NbFkmSqsgB3co3b6tMgg/dCLedA7PPySt0j54A2x6Zb9xWXdtBpaqdQlt+veEDsOvJyx9rpK5QZSnBKmtVlnbsxMEtiyRJVWRwoazQlgds7/pB2O1Dyx9LRQMLDY5ie772iu25i9DwMXn61fL+rrq7DhulK1RnhbY81umy/+i7a9QOxwxNmSRJqgKDCy2vu5u3VggsugugUhEIu4INllTMLWQ3/QjuvBAWvwxtI2DrQ2GPf4d1+xgH1OhGjM2B/F9O6znNqAmw2yl5tigHdUuSGoDBhZSK8NI8+Pv/wV0XwcL5MH492PFYmHoCjFilMZ+OdycVO2b8ikLtAsdiEebdAeccnmdEKmtfAv+8CO65FN45HbapcCazRpQS7PdFeO2F7qejHTsR3vNrWHVSawT4kqSm0FLBRUQMBz4C7AjsBGwNDAdOTilN72dek4FHe0lyYUrp6JUrqYZMKsJjf4NfvTs/OS977kG4+ivwj+lw/B9h/AaNHWCUuxm9/C949Pr8fv1dYeIWPXdBGlQJLnzv8oFFZ+1L4eIP5pnMxqzZnK1HEZCAt38Pdv9wxwrdw0blqXe3OwoKrtAtSWosLRVcAGOBH5R+fhr4F7DBAPO8A7ikm/13DTBfDbaUckBxwTHLBxadLXgCzn83fPSWoS1bNRXbYeFz8IdPwANXLL9WyUZ7wtu+n4OMzjexnQOOao+5KbbDvX+Al57qPd2yxXDrGbD3f1Tvd9ebctC0+iZw4Lc69tck4JMkaeBaLbhYCBwM3J5SmhcRXwW+MsA8b08pfXWgBVMtJJh9Nixa0HuyZ++Dh6+FjfdpvBu+YjGf3y/fAi/MWfH4YzfBGW+Fk67ON7hRyDe291ySz7m4FNbZAXZ6X56KuBoKbfDQVZWlffAq2Oez1fm99azrddVo15kkSSUt1d6eUlqSUro8pTSv1mVRHYhCfoJeiXt+35g3fIUC3PDd7gOLskUvwpVfLM3aVIS2YXktidfmw52/hj9/EU7bIo9JqZb2JRWmW1q93ylJkgZdSwUXg2TdiDglIr5Q2vawjLDq0pJXK0u3tMJ09WbZErj9/L7TPXgVvPg4PHw1XPJhSO1w9Pmw8/tL+SyCy/5fDjY6d6taGSn1vNp2V5O2HdjvkiRJQ8rgYuAOAH4OfKO0vSMirouIDWtbLFVkzc0rS7dGhenqzUtP9TxourNUhGfuhVXWzsHIL/aHh66Gt38f1tulI93Mb644NqO4LAcM/bHT+/LA5b684eTGXCRPkqQWZXCx8hYCXwd2AVYrvfYBrgOmAddExNhKMoqIWd29gCmDU3QB+aZ1lxP6TlcYBru8v/830PWgbXg/0o7IYywgt1T85oQ80H33j3SkWfxSnjoVYOlCePgauP/yjm5XlQQCEXn19/2/3Hu6nd+f17poxO5okiS1qIYLLiJiTkSkfrzOHYxypJSeSSl9OaU0O6X0Yul1A/AW4BZgM+CkwfjdqpJCG2yyD2zx1t7T7fFv+Yl+I06HOn59mLhl3+lGjYcNdoWn7+7Yt+QVuP28vKjdmDVg9GpwwuW5xeHyz+ZxGOcdBRceCz/aEc46NLd+VNptao+PwqE/ztP8djZmddj3C3DIDwbeBUuSJA2pRpwt6mFgUT/Szx2sgnQnpbQsIqYDuwF7Az+s4DO7dLe/1Hqxc3VLqOUU2+Fd58KfPgV3XrD8AOIRq8Ce/wbTPp9bLRoxuEhF2PWD8KdP955ux/fmFaNvPWP5/Q9elYOANbeAnY6FCRvBOe+Ax/66Yh6PXp9nnjr+TzBpu8paHHY6Nv/uR2fmLlyjV4fNDoBhIxq3ziVJamENF1yklPavdRkq8GxpW1G3KNVQoQ0IOOzH8OavwF2/La3QvQFs+858w93IN7lRgKknwuM3wz9/032ajfaE/U+FR2/Iq2Z3Vp7VacwasN2Reere7gKLsiWvwJ8+CSdfV2H5bu5maAAAFxRJREFUIr823oe8ohy5G1r5mCRJaigNF1w0iN1L20dqWgpVplDqHThmzfyUv7zWQ/nJezPc5B4xHTbbP08n+9TsvG/NzWHqB2DqCXmmqIu6GX8yabuO7bBReRXpvjw1G+bdmaezrXS8hOMqJElqCgYXfYiI8cA6wILO62NExG7AbSmlJV3S7wd8svR2UMZ7aJBEAKVAopludsuzO233LtjhmDwQu1iEkavklom7fgdXfK5joHZnU0+EJ2+F4WPyitnP3FvZ73xqFqy9dfXOQZIkNYSWCy4i4nN0zMK0Y2l7QkS8sfTzjSml6Z0+8g7gTOAs4PhO+78NbBMRM4EnS/u2B/Yr/XxqSumm6pZeGoBywDR8TN4+dDVcfAq8+tzy6dbaGt7wgRyIjBibA5F1dli53yVJklpKywUXwIHkKWM727P0KptO384hBx5vAA4ChgNPA78GfpxS+svAiyoNos3eDHv/P7jx+/Dyv3LLzZu/Bnt9PLdS3H0xPP9QDjB2eA+sujZstFfvYy4gt5RsMo0GnIxu5XXuRtfbPkmSmlzLBRcppWn9TD8DmNHN/l8CFXRAl+rYbh+CN5wED10D49bLK2LfegZc85/Ld5O6/jvwmQfymJS+govND4AJLbKGZHmw/9zb8niWJ/6e32+0F+x2SseYFUmSWkTLBReSuigMy60YhbYcWPzxkyumWbow3zy/6dPwyPEwa0b3ea02GQ49vXWe2kfk8So3/2z5/fMfgdvOgWmfy1MZS5LUIlqo34KkHhXaYNmS3GLRk+v+O697ccgP4bCfLD8OY/RqsOe/5ylox6zZGoFFKsI/pq8YWHQ281vwz4sqW7lckqQmYMuFpHzze88l3c8Y9XqaZXDBMXDYT/Ng752OzWM1iu2wylrQNjzfcEeLPLOIAvztJ32nu+n0vEaIJEktoEXuAiT1qtCWB2/3pX1p7jZVaIN7/wgPXglLX+u0Jsgg/kkpFru8r3FrwNzbc/envsy7Pa8jIklSC7DlQlJucRixSmVpR5Smsn30+jwO46DvwOobv75EyKCULQrwwiMw6yxY8ASMGg/bHgEb712b8R0pwaIFlafvT1pJkhqYwYWkbOtD4eov5xvn3mx1aN4+eWve9ncNjP4oFiG1w+8/DHdcsPyxWTNgnR3hPRfC2IlDG2BE5MHrlSgMg/HrD2pxJEmqF3aLkpRbBlabDFsc2Hu6tuGw68nw1GyYOzsvuLfh7oN3Y18o5G5YXQOLsnm3w9mH1qaL1Gob5ZaTvmx5cB7wLklSCzC4kJQV2+Hwn/fcEtE2PB+fOAX+8j8wbCQcfFrutjQYUhEWPAm3n9d7umfvh7tqMCNTsR32/WKul54MHwPTPlv78SGSJA0RgwtJWaENRq4KH7gKDvkRrLszjBwHq07KC+196MY869FVX4aX5sJxv4fJew3eIO4o5MCikuDl9vOGftxFoS232hx9QZ4tq6vx68Oxv4O1tmmNqXklScIxF5I6K7Tl107Hwi7vX/7Yq8/DozfAtkfCAf85NE/jX5pb3XSDYbP94VP3wj2/hyduyeMxJr8pdzGLyC9JklqEwYWkFXV90p4SjB4PG+3J6w2eQ/E0fswaFaZbfXDL0ZsIiGGw9WGw7TvzvlZZoVySpC7sFiWpbxF51qPCsDzIeigU22H7d1eWdrt3Dd7Yj0p1DiYMLCRJLcrgQlJ9KrTBxC37nsFq1Um5G1errAwuSVId89tYUv0qtsORZ8DkN3Z/fNy68L5LYdiooS2XJEnqlmMuJNWvQlsOHI7/Ezx0Dcw+G158vGOF7u3fVeqqZTckSZLqgcGFpPpWDhw2mZZnZipz0LQkSXXHblGSGkPXQMLAQpKkumNwIUmSJKkqDC4kNZbOi/cNxUJ+kiSpYo65kNQYymMs5t0B918G7Utgra1hm3fCsBG1Lp0kScLgQlIjKLbDS0/Bb46Hp2Ytf+yKz8FbvwE7vrcmRZMkSR0MLiTVt5Rg0QI48yBY8OSKx197AS75CLSNgG2PzKuJS5KkmnDMhaT6FgG3/Lz7wKKza/4TSENSJEmS1D2DC0n1b/bZfad58XF4+DoHeUuSVEMGF5Lq27LF8PK8ytLOfwRbLyRJqh2DC0n1rW14Hk9RiZGrDG5ZJElSrwwuJNW3KMBWh/Sdbvho2PJgCFfuliSpVgwuJNW3Yjvs/pEcZPRmh/fAqPHOFiVJUg0ZXEiqb4U2WH8qHPLD/HN3Nt0fDvwWFItDWzZJkrQc17mQ1Bh2Pg423AP+/r/wwJ9h2SJYaxuYeiJMeTuQoODzEkmSasngQlLjWH0TOPi0/CorthtUSJJUJ/xGltQ4uusW1VNXKUmSNOQMLiRJkiRVhcGFJEmSpKpoqeAiIjaPiM9GxLUR8URELImIpyPi0ojYdyXz3DMiLouI+RGxMCLujIhPRDjZviRJklpLSwUXwNeBbwFrA5cB/wP8FXgbcG1EfKw/mUXEYcANwN7AxcBPgBHA94ELqldsSZIkqf612mxRVwDfTind1nlnROwDXAV8NyJ+k1Ka11dGETEO+AXQDkxLKd1a2n8qcC1wZEQcnVIyyJAkSVJLaKmWi5TSjK6BRWn/9cBMcqvDnhVmdyQwEbigHFiU8loEfKn09sMDKrAkSZLUQFoquOjD0tJ2WYXp9yttr+jm2A3AQmDPiBg50IJJkiRJjcDgAoiIjYD9yQHBDRV+bMvS9oGuB1JKy4BHyd3ONqlGGSVJkqR612pjLlZQalk4DxgJ/EdK6YUKPzq+tF3Qw/Hy/gkVlGFWD4emVFgWSZIkqeYaruUiIuZEROrH69xe8moDzgH2Ai4ETqtmUUvbVMU8JUmSpLrViC0XDwOL+pF+bnc7S4HFucBRwK+BY1NK/QkEyi0T43s4Pq5Luh6llHbpoYyzgJ37USZJkiSpZhouuEgp7T/QPCJiGHA+ObA4HzgupdTez2zuB6YCWwDLdWsq5b8xeXD4IwMtryRJktQIGi64GKiIGEFuqTgMOBs4IaVUXImsrgXeCxwI/KrLsb2BMcANKaXFAyju5HvvvZdddum2YUOSJEmqinvvvRdg8kDzif71BGpspcHbvwMOBn4JfLCvwCIixgPrAAs6L65XWkTvYXL3p706LaI3ihx47AEcM5BF9CLi0VL+c1Y2jyopDyy/r6alaAzWVWWsp8pZV5WxnipnXVXOuqqM9VS5eq6rycBLKaWNB5JJqwUXZwLHA88BP6X7wdYzU0ozO33meOBM4KyU0vFd8jscuIg8BuQCYD5wKHma2ouAd/VzHEddKs9m1dPYEHWwripjPVXOuqqM9VQ566py1lVlrKfKtUJdtVq3qHIktibw5V7Szawks5TSJRGxD/BF4AhgFPAQ8CngR80QWEiSJEmVaqngIqU0bSU+MwOY0cvxv5K7WUmSJEktreHWuZAkSZJUnwwuJEmSJFWFwYUkSZKkqmip2aIkSZIkDR5bLiRJkiRVhcGFJEmSpKowuJAkSZJUFQYXkiRJkqrC4EKSJElSVRhcSJIkSaoKgwtJkiRJVWFwoddFxOYR8dmIuDYinoiIJRHxdERcGhH7rmSee0bEZRExPyIWRsSdEfGJiGirdvmHSkQMj4iPR8SZEXF7qZ5SRJy0EnlNLn22p9cFg3EOQ6WaddUpz6a7pjqr1vn1cV3dPFjlr6aIWD8izoiIuRGxOCLmRMQPImK1fuazeulzc0r5zC3lu/5glX0oVaOeImJmH9fMqME8h6EQEUdGxOkR8ZeIeKl0XueuZF5VuTbrUbXqqVQnPV1P/xqMsg+liFgjIk6KiIsj4qGIeC0iFkTEjRHxgYjo1z12M11Tw2pdANWVrwPvBu4BLgPmA1sChwKHRsTHU0o/qjSziDgM+C2wCLiwlN8hwPeBvYCjqlr6oTMW+EHp56eBfwEbDDDPO4BLutl/1wDzrbWq1lUTX1PAoJzfY8CMbvY/ufKlHBoRsSlwE7AWcClwH7Ar8HHgwIjYK6X0fAX5rFHKZwvgWuACYApwAvC2iNgjpfTI4JzF4KtWPXXytR72LxtQQevDl4AdgFfI/wemrEwmg1Dn9aYq9VSygI7vgM5eGUCe9eIo4GfAPOA64HFgbeCdwHTgoIg4KlWwWnXTXVMpJV++SCkBHA/s1M3+fYAlwGJgnQrzGgc8U/rM1E77R5H/AyXg6Fqf80rW0wjgoHJdAF8tnc9JK5HX5NJnZ9T6vBqgrpr2mhqM8yuln1nr8xpAffy5dA7/3mX/90r7f15hPv9bSv+9Lvs/Vtp/Ra3PtU7qaWa+Jaj9OQ1iXe0LbA4EMK1UP+fWqs7r9VXFepoDzKn1+QxiPe1HfvhT6LJ/EjnQSMARrXhN2S1Kr0spzUgp3dbN/uvJXzwjgD0rzO5IYCJwQUrp1k55LSI/FQH48IAKXCMppSUppctTSvNqXZZ6V+W6atprqqTZz69iEbEJ8BbyzclPuhz+CvAq8L6IGNtHPmOB95XSf6XL4R+X8n9r6fc1nGrVU6tIKV2XUnowle7aVkYr1Hk16qkVpJSuTSn9IaVU7LL/X8DPS2+n9ZVPM15TBheq1NLSttKm8f1K2yu6OXYDsBDYMyJGDrRgTWLdiDglIr5Q2m5f6wLVoWa/pgbj/CZExIml6+qjEbH7gEs5NMp1cWU3X9wvA38FxgB9nc8ewGjgr6XPdc6nCFxZertSY8rqQLXq6XUR8e6I+FxEfCoiDmrg/0+Dpep13uRGRsSxpb9BH4+Iffs7fqxB9eeeqemuKcdcqE8RsRGwP/nm5oYKP7ZlaftA1wMppWUR8SiwDbAJcG81ytngDii9XhcRM4H3p5Qer0mJ6k+zX1ODcX47AL/svCMi7gDel1L65wDKOth6rIuSB8lP+rYArhlgPpTyaUTVqqfOuk4i8UxEfDSldNFKlK8ZDUadN7NJwDld9j0aESeUekU0nYgYBhxXetvdw6Kumu6asuVCvSo9tToPGAl8NaX0QoUfHV/aLujheHn/hAEUrxksJA+k3wVYrfTahzw4bBpwTSM1hQ6yZr+mqn1+3yMPAp8IrAq8AbiIHHBcGxHrrWQ5h0K16sJrJqvk/C4l9x9fn9zaMwX4ZumzF0bEQQMoZzNp9muqms4kP5icRJ7cYzvyGKjJwOURsUPtijaovgVsC1yWUvpzBemb7poyuGgyfUz91t2rx+nlSk2X55BvUC4ETqtmUUvbmvTprGY9DURK6ZmU0pdTSrNTSi+WXjeQn1LcAmwGrPS0rdVQL3VVSVFL25r1Ex7kuurX+aWUPp1Suiml9FxK6ZWU0q0ppaPIs1GtCXymn6dXT6r1b13za2aQVXx+KaXvp5T+mFJ6KqW0KKV0f0rpC8CnyfcK/z2YBW0izX5NVSyl9LXSuISnU0oLU0p3pZQ+RH7wMZo8wUdTiYiPkf/P3Ece71WVbEvbhrmm7BbVfB4mT2NZqbnd7SwFFueSp1r7NXBsPwd3lSPt8T0cH9cl3VCrSj0NllI3mOnAbsDewA+H8vd3US91Ve/XFAysrobq/H4OHEG+rupVteqiEa6ZgRiK85tOngp5x4hYtevYlRbU7NfUUPg5+Qa8nv8G9VtEfJT8XX0PsH9KaX6FH226a8rgosmklPYfaB6l/oLnkwOL84HjUkrt/czmfmAquY/grG7y35g80Kkm88tXo56GwLOlbU27RdVRXdX1NQUDrquhOr+6uK76cH9p29NYiM1L2576KFc7n3o16OeXUloUES+Tu2yOBVo9uGj2a2ooPFPa1vPfoH6JiE+Qg/C7yIHFM318pLOmu6bsFqXlRMQIcr/so4CzyQM/+xtYQF6sCuDAbo7tTZ754KaU0uKVKmhrKM8M0bALfFVZs19TQ3V+jXBdXVfaviW6rHIbEauSu2q+BvS10vjNpXR7lT7XOZ8Cufth59/XaKpVTz2KiC3JgcXLwHMrm08TGfQ6bwF7lLb1/DeoYhHxWXJgcTuwbz8DC2jCa8rgQq8rDd6+GDiMPMPMCV2nRevmM+MjYkpErNPl0EXkL6KjI2Jqp/SjgP8qvf1Z1Qpf53qqp4jYrRTQdU2/H/DJ0ttajWGoiRa+pvp9fhExplRXG3bZv3N3EwFEnuL4G6W3dXtdpZQeJk8TOxn4aJfDXyM/8Tw7pfRqeWepHpZbSTil9Ap53NhYVuzf/W+l/P+cGnSF7mrVU0Rs0t0A/4hYkzwoF/L6K82wSndFImJ4qa427bx/Zeq8mfVUTxGxTUSs3k36jchrzEAd/w2qVEScSh7APYvcYtFjAN5K11S4RorKIuJM8irdzwE/pfvBQzNTSjM7feZ48pfPWSml47vkdzj5hmkReXrD+cCh5GnXLgLe1aiL9ETE58izqQDsSJ6B5yY6pra8MaU0vVP64+mmniJPN7sNeZHCJ0u7t6dj3utTU0rlG8uGVK26Kh1r2msK+n9+ETGN/NTr+pTStE77ZwDvJLeGPEFe9XsKuVWkDfgFcEo911XpC/gmYC3yTEb3kscg7UvuHrBnSun5TunzsuQpRZd81ijlswW5Pv4ObEV+iPJMKZ+HB/t8Bks16qn0f246cD153NB8YEPgYHI/8FuBA1JKLw7+GQ2e0v+vw0tvJwFvJT89/0tp33Mppc+U0k4GHgUeSylN7pJPv+q80VSjniLiq8DnyH+fHiW3fG0KvA0YBVwGvCOltGRQT2YQRcT7gRlAO3A63Y+JmJNSmlFKP5lWuaZSHSwT7qs+XuQb3NTH66tdPnN8af+MHvLci/xH5AVys94/yU/k22p9voNcVzO6pO+2noAPAH8kr8z5Cvkm8HHy7FxvqvV51lNdNfs1tTLnR56uOJGD/s77Dwd+BzwEvAQsAeYBfwAOrfU59qMuNiAHmvNK5/AYecDk6t2kTfkrrdt8Vi997rFOdXEGsH6tz7Ee6ok8ReiM0rX2PHkBsPnkm8l/B0bU+hyrVE9f7eNv0ZxOaSd33beydd5or2rUE3lK9V+RZ016sXRNPQtcRV4DImp9nkNQT8v9bW6la8qWC0mSJElV4ZgLSZIkSVVhcCFJkiSpKgwuJEmSJFWFwYUkSZKkqjC4kCRJklQVBheSJEmSqsLgQpIkSVJVGFxIkiRJqgqDC0mSJElVYXAhSZIkqSoMLiRJkiRVhcGFJEmSpKowuJAkNaWI+EBE/G9E3BIRCyMiRcR/1bpcktTMhtW6AJIkDZL/AcYDLwBzgU1rWxxJan62XEiSmtXRwOSU0uqALRaSNAQMLiRJDSEiLil1bfr3bo59vXRsenlfSumKlNJjQ1tKSWptBheSpEZxIvA48N2I2Km8MyL2B74A3AN8rEZlkyRhcCFJahAppfnAMUAbcGFErBIRawHnAouBd6WUFtayjJLU6gwuJEkNI6V0E3AqsDnwv+TAYhLwsZTS3bUsmyTJ2aIkSY3n28A04D2l979KKU3vObkkaajYciFJaigppQRc3GnXD2pVFknS8gwuJEkNJSI2B04jr19RBKZHxKjalkqSBAYXkqQGEhEjgQuBseR1LL4JbIetF5JUFwwuJEmN5DRgJ+A7KaUrga8AfwVOiYh31bRkkiQid12VJKm+RcTh5LEWtwBvTCktK+3fALidPEnJTimlR0r7TwLeWPr4ZsBewJ3AbaV996WUvjV0ZyBJzc/gQpJU9yJiQ3IAUSAHEI92OX4YcAnwD3LgsSQiZgDv7yXb61NK0wanxJLUmgwuJEmSJFWFYy4kSZIkVYXBhSRJkqSqMLiQJEmSVBUGF5IkSZKqwuBCkiRJUlUYXEiSJEmqCoMLSZIkSVVhcCFJkiSpKgwuJEmSJFWFwYUkSZKkqjC4kCRJklQVBheSJEmSqsLgQpIkSVJVGFxIkiRJqgqDC0mSJElVYXAhSZIkqSoMLiRJkiRVxf8HvR3EOhiGUcgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "image/png": {
       "height": 263,
       "width": 395
      },
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "seed(314)\n",
    "nobs = 30\n",
    "x1 = multivariate_normal([-1, -1], np.eye(2) * 0.3, nobs)\n",
    "t1 = np.ones(nobs)\n",
    "\n",
    "x2 = multivariate_normal([1, 1], np.eye(2) * 0.3, nobs)\n",
    "t2 = np.zeros(nobs)\n",
    "\n",
    "D = np.c_[np.r_[x1, x2], np.r_[t1, t2]]\n",
    "D = pd.DataFrame(D, columns=[\"x1\", \"x2\", \"t\"])\n",
    "\n",
    "sns.scatterplot(\"x1\", \"x2\", hue=\"t\", data=D)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0., 1.],\n",
       "       [0., 1.],\n",
       "       [0., 1.],\n",
       "       [0., 1.],\n",
       "       [0., 1.]])"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ohe = OneHotEncoder(categories=\"auto\", sparse=False)\n",
    "T = ohe.fit_transform(D[\"t\"].values[:, np.newaxis])\n",
    "T[:5]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1.    , -0.909 , -0.5717],\n",
       "       [ 1.    , -0.5332, -1.3873],\n",
       "       [ 1.    , -1.5103, -0.5144],\n",
       "       [ 1.    , -1.1215, -0.7909],\n",
       "       [ 1.    , -1.4232, -0.5273]])"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Xtilde = np.c_[np.ones(2 * nobs), D[[\"x1\", \"x2\"]]]\n",
    "Xtilde[:5]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(3, 2)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "W = inv(Xtilde.T @ Xtilde) @ Xtilde.T @ T\n",
    "W.shape # (nfeats, nclasses)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Up to this point, it is easy to see that the constraint\n",
    "\n",
    "$$\n",
    "    {\\bf a}^T {\\bf t}_n + b = 0\n",
    "$$\n",
    "\n",
    "is satisfied for a one-hot-encoded vector ${\\bf t}_n$ by choosing $\\bf{a} = \\bf{1}$ (a vector of ones) and $b = -1$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.0"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a = np.ones(2)\n",
    "b = -1\n",
    "a @ T[0] + b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Or more generally,\n",
    "np.all(T @ a - 1 == 0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, considering a new set of points generated by a different probability distribution,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1.    ,  0.0042, -0.2876],\n",
       "       [ 1.    ,  1.0095, -0.2514],\n",
       "       [ 1.    ,  0.3646,  1.1981],\n",
       "       [ 1.    ,  0.2661,  0.3501],\n",
       "       [ 1.    ,  0.2703, -0.326 ]])"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Xtest = multivariate_normal([0, 0], np.eye(2) * 0.3, nobs)\n",
    "Xtest = np.c_[np.ones(nobs), Xtest]\n",
    "Xtest[:5]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.4452, 0.5548],\n",
       "       [0.7067, 0.2933],\n",
       "       [0.8043, 0.1957],\n",
       "       [0.6265, 0.3735],\n",
       "       [0.5058, 0.4942]])"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# We can make batch estimations via\n",
    "yest = Xtest @ W \n",
    "yest[:5]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,\n",
       "        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0., -0.,\n",
       "        0., -0.,  0.,  0.])"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Which, as we see, does satisfy the imposed linear constraint\n",
    "yest @ a + b"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To understand why this is the case,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 87,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1.]])"
      ]
     },
     "execution_count": 87,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "-b * np.ones((1, 2 * nobs)) @ Xtilde @ inv(Xtilde.T @ Xtilde) @ Xtest[[0]].T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 105,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1., -0.,  0.]])"
      ]
     },
     "execution_count": 105,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.ones((1, 2 * nobs)) @ Xtilde @ inv(Xtilde.T @ Xtilde).T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 102,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.0166, -0.0163,  0.0051],\n",
       "       [ 0.0173,  0.0172, -0.0289],\n",
       "       [ 0.0162, -0.0383,  0.022 ],\n",
       "       [ 0.0166, -0.018 ,  0.0038],\n",
       "       [ 0.0163, -0.035 ,  0.0194],\n",
       "       [ 0.017 , -0.001 , -0.017 ],\n",
       "       [ 0.0164, -0.0285,  0.014 ],\n",
       "       [ 0.0171,  0.0036, -0.0184],\n",
       "       [ 0.0167, -0.0013,  0.0001],\n",
       "       [ 0.0167, -0.0018, -0.0034],\n",
       "       [ 0.0165, -0.0349,  0.0098],\n",
       "       [ 0.0157, -0.0597,  0.0489],\n",
       "       [ 0.0168, -0.0079, -0.0048],\n",
       "       [ 0.017 ,  0.0102, -0.0165],\n",
       "       [ 0.0174,  0.0162, -0.0334],\n",
       "       [ 0.017 , -0.0016, -0.0167],\n",
       "       [ 0.0171,  0.0029, -0.0212],\n",
       "       [ 0.0162, -0.0349,  0.0223],\n",
       "       [ 0.0165, -0.0275,  0.009 ],\n",
       "       [ 0.0168, -0.0172, -0.0047],\n",
       "       [ 0.0164, -0.0307,  0.0133],\n",
       "       [ 0.0171,  0.0042, -0.0203],\n",
       "       [ 0.017 ,  0.0059, -0.0143],\n",
       "       [ 0.0174,  0.0288, -0.0353],\n",
       "       [ 0.0171,  0.002 , -0.0205],\n",
       "       [ 0.017 , -0.0038, -0.0148],\n",
       "       [ 0.0166, -0.0202,  0.0044],\n",
       "       [ 0.0174,  0.0197, -0.0336],\n",
       "       [ 0.0171,  0.0022, -0.0198],\n",
       "       [ 0.0173,  0.012 , -0.0285],\n",
       "       [ 0.0162, -0.004 ,  0.0237],\n",
       "       [ 0.0163, -0.0014,  0.015 ],\n",
       "       [ 0.0163, -0.0041,  0.0158],\n",
       "       [ 0.0162, -0.0093,  0.0236],\n",
       "       [ 0.0157, -0.0363,  0.0469],\n",
       "       [ 0.016 , -0.0153,  0.0301],\n",
       "       [ 0.0162, -0.0159,  0.022 ],\n",
       "       [ 0.0167,  0.0112, -0.0034],\n",
       "       [ 0.0169,  0.0316, -0.0109],\n",
       "       [ 0.0157, -0.0322,  0.047 ],\n",
       "       [ 0.017 ,  0.042 , -0.019 ],\n",
       "       [ 0.0166,  0.0127,  0.001 ],\n",
       "       [ 0.0165,  0.002 ,  0.006 ],\n",
       "       [ 0.0176,  0.0565, -0.0443],\n",
       "       [ 0.0171,  0.0336, -0.0216],\n",
       "       [ 0.0173,  0.0441, -0.0295],\n",
       "       [ 0.0168,  0.0218, -0.0085],\n",
       "       [ 0.0164,  0.0039,  0.0131],\n",
       "       [ 0.0168,  0.0275, -0.0083],\n",
       "       [ 0.0161, -0.0162,  0.0259],\n",
       "       [ 0.0171,  0.0343, -0.0236],\n",
       "       [ 0.0165,  0.0078,  0.0056],\n",
       "       [ 0.0165,  0.0132,  0.0057],\n",
       "       [ 0.0159, -0.0181,  0.0384],\n",
       "       [ 0.0164,  0.0023,  0.0133],\n",
       "       [ 0.0166,  0.003 ,  0.0024],\n",
       "       [ 0.0174,  0.0577, -0.0365],\n",
       "       [ 0.0159, -0.0188,  0.0362],\n",
       "       [ 0.0166,  0.0201,  0.0026],\n",
       "       [ 0.0164, -0.0001,  0.0114]])"
      ]
     },
     "execution_count": 102,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Xtilde @ inv(Xtilde.T @ Xtilde)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 99,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.0174],\n",
       "       [ 0.0415],\n",
       "       [-0.0337]])"
      ]
     },
     "execution_count": 99,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "inv(Xtilde.T @ Xtilde) @ Xtest[[1]].T"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
